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SECTION  1 
INTRODUCTION 


The  Protocol  to  the  Threshold  Test  Ban  Treaty  (TTBT)  provides  for  verifica¬ 
tion  of  a  150  kt  limit  on  the  yield  of  underground  nuclear  tests  (UGT).  The  Protocol 
provisions  for  on-site  inspection  (OSI)  include  yield  estimation  by  active  close-in 
measurements  during  an  UGT.  The  Defense  Nuclear  Agency  (DNA)  has  the  respon¬ 
sibility  for  yield  estimation  on  a  class  of  tests  defined  as  'non-standard  in  the 
Protocol,  and  has  developed  a  method  called  'HYDROPLUS'  for  analyzing  these 
tests.  HYDROPLUS  provides  a  yield  estimate  by  comparing  measured  values  of  the 
shock  stress,  particle  velocity,  and  time-of-arrival  (or  shock  velocity)  in  the  geologic 
medium  surrounding  the  test  region  with  theoretical  predictions  of  these  quantities. 

In  the  HYDROPLUS  methodology,  the  theoretical  predictions  are  based  on 
first-principles  calculations  of  the  evolution  of  the  explosion-produced  shock  wave 
in  the  surrounding  medium.  Accurate  predictions  require  the  use  of  an  appropriate 
computer  code  to  solve  the  equations  of  motion,  an  adequate  idealization  of  the  test 
configuration,  and  accurate  material  models  for  the  geological  media  and  other 
materials  surrounding  the  source. 

For  Treaty  verification,  yield  uncertainties  better  than  30%  with  high  confi¬ 
dence  are  required.  This  means  that  the  combined  error  due  to  measurements  and 
predictions  must  not  exceed  30%.  This  places  stringent  requirements  on  the  accuracy 
of  both  the  measurements  and  the  predictions. 

The  development  of  accurate  equations  of  state  (EOS)  for  geologic  materials 
has  emerged  as  a  critical  accomplishment  for  the  success  of  HYDROPLUS.  Early  tests 
of  the  methodology  revealed  that  conventional  equilibrium  equations  of  state  were 
not  capable  of  accurate  predictions  of  the  stress,  particle  velocity,  and  shock  velocity 
simultaneously.  Consistent  predictions  of  all  three  shock  features  was  shown  to 
require  non-equilibrium  modeling  of  the  phase  transformation  between  low 
pressure  and  high  pressure  rock  phases. 

The  application  of  HYDROPLUS  also  required  the  development  of  accurate 
site-specific  equations  of  state  on  a  short  time  line  consistent  with  Protocol  require¬ 
ments. 

S-Cubed  has  developed  a  parametric  computational  equation  of  state  that 
accurately  represents  the  observed  non-equilibrium  behavior.  This  equation  of  state 
has  been  implemented  as  part  of  a  computer  workstation  encapsulation  called 
HYPAWS,  for  HYDROPLUS  Analysis  Workstation.  This  program  provides  a 
convenient  interactive  capability  for  comparing  the  theoretical  equation  of  state 
with  laboratory  test  data,  and  varying  model  parameters  to  achieve  the  best 
agreement  between  the  theoretical  model  and  the  data. 
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This  report  reviews  the  theoretical  foundation  for  the  non-equilibrium 
equation  of  state  for  silicate  rocks,  and  documents  the  parametric  computational 
equation  of  state  developed  by  S-Cubed  for  the  DNA  HYDROPLUS  methodology. 
Section  2  presents  a  qualitative  discussion  of  the  requirement  for  non-equilibrium 
equations-of-state,  and  the  role  of  theoretical  modeling  in  the  formulation  of  the 
computational  equation  of  state.  Section  3  presents  a  qualitative  discussion  of  the 
non-equilibrium  model,  while  section  4  provides  the  mathematical  details  of  the 
model.  Section  5  presents  additional  details  on  the  code  implementation  of  the 
model  and  illustrates  the  results  with  an  application  to  rhyolite.  Conclusions  and 
recommendations  are  presented  in  Section  6. 
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SECTION  2 

EOS  MODELING  BACKGROUND 


2.1  NEED  FOR  A  NON-EQUILIBRIUM  EQUATION  OF  STATE. 

For  some  time,  laboratory  shock  wave  data  have  shown  that  granite  and 
rhyolite  display  release  paths  that  fall  below  the  Hugoniot  on  a  P-V  plot.  Such 
release  patfe  are  known  to  be  inconsistent  with  equilibrium  thermodynamic 
behavior. 

Release  paths  have  a  direct  influence  on  shock  wave  attenuation,  and  shock 
wave  measurements  on  early  tests  of  the  HYDROPLUS  method  confirmed  that 
releases  beneath  the  Hugoniot  occurred  in  full  scale  tests  as  well  as  laboratory  scale 
tests.  Consistent  with  this,  it  was  demonstrated  that  equilibrium  equations  of  state 
were  not  able  to  reproduce  all  the  shock  wave  data. 

The  observed  non-equilibrium  behavior  has  a  natural  physical  explanation 
associated  with  the  non-equilibrium  composition  of  polymorphic  (distinct  solid) 
phases  of  silica.  This  non-equilibrium  phase  phenomenon  is  characterized  as  phase 
hysteresis  and  profoundly  affects  groimd  shock  propagation.  The  presence  of  phase 
hysteresis  in  large  scale  explosions  has  only  become  appreciated  in  the  last  few  years. 
Previously,  the  best  computational  EOS  models  assumed  phase  equilibrium 
[Andrews,  1971;  Lee  et  al.,  1984]. 

The  inadequacy  of  equilibrium  EOS  models  emerged  during  the  Defense 
Nuclear  Agency's  HYDROPLUS  program.  The  HYDROPLUS  requirement  to  make 
accurate  predictions  of  particle  velocity  and  shock  pressure  in  addition  to  time-of- 
arrival  made  imprecedented  requirements  on  the  accuracy  of  computational 
models.  The  failure  of  conventional  models  to  accurately  predict  all  measured  shock 
properties  was  apparent  to  all  participants  in  this  program.  Schuster  [1993],  using 
empirical  models,  was  the  first  to  note  that  hysteretic  models  that  agreed  with 
laboratory  data  also  produced  groimd  shock  predictions  which  agreed  well  with 
HYDROPLUS  data.  This  demonstrated  that  hysteretic  unloads  were  not  just  an 
artifact  of  laboratory  experiments,  but  were  essential  to  get  detailed  agreement  with 
ground  shock  data. 

The  first  ground  shock  calculations  based  on  a  thermodynamic  model  of  non¬ 
equilibrium  phase  transitions  were  reported  by  Swegle  [1990].  This  paper  contains  a 
useful  review  of  the  evidence  for  non-equilibrium  shock  behavior  of  silica,  the 
formulation  of  a  theoretically  sound  computational  model,  and  calculations  which 
demonstrate  the  sensitivity  of  time-of-arrival  to  model  details.  Similar  non¬ 
equilibrium  models  were  developed  independently  by  S-Cubed  [Baker  et  al.,  1993] 
and  other  DNA  HYDROPLUS  calculators. 
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2.2  THEORETICAL  EQUAHONS-OF-STATE. 


While  purely  empirical  (based  on  data  fitting  with  no  theoretical  constraints) 
EOS  models  have  had  significant  successes  in  ground  motion  predictions,  they  can 
be  used  with  confidence  only  for  a  limited  range  of  conditions.  Important  physical 
processes  such  as  melting,  vaporization,  chemical  dissociation,  and  ionization 
cannot  be  accurately  modeled  by  purely  empirical  methods.  These  processes  must  be 
realistically  treated  in  the  EOS  for  proper  behavior  at  high  temperatures  and  pres¬ 
sures.  Because  of  this,  computational  equations  of  state  (EOS  subroutines)  should 
have  a  good  theoretical  foundation. 

Theoretical  EOS  models  of  the  various  physical  processes  have  been 
implemented  independently  by  a  number  of  agencies.  The  documentation  of  the 
PANDA  n  code  [Kerley,  1991]  probably  provides  the  best  overview  of  these  methods. 
The  thermal  component  of  the  HYPAWS  EOS  is  based  on  analogous  models 
developed  at  S-Cubed. 

The  most  significant  element  of  the  HYPAWS  equation  of  state,  however,  is 
the  effect  of  non-equilibrium  solid-solid  polymorphic  phase  transitions.  Again, 
theoretical  considerations  provide  improved  confidence  over  purely  empirical 
methods.  In  particular,  thermodynamics  is  used  to  assure  physically  acceptable 
behavior.  An  important  consequence  of  the  thermodynamic  models  is  that  the 
hysteretic  unloading  behavior  observed  in  laboratory  experiments  follow  naturally 
from  the  second  law  of  thermodynamics,  without  curve  fitting. 

The  HYPAWS  EOS  integrates  a  non-equilibrium  thermodynamic  model  of 
polymorphic  phase  transitions  and  an  (equilibrium)  thermal  model  of  high 
temperature  processes  in  a  consistent  manner.  However,  in  contrast  to  the  high 
temperature  thermal  model,  the  non-equilibrium  phase  model  needs  empirical  data 
to  calibrate  its  model  parameters.  It  bears  emphasis  that  the  resulting  model  satisfies 
the  general  theoretical  requirements  of  thermodynamics,  whereas  a  purely 
empirical  model  will  not. 

This  semi-empirical  approach  is  necessary  for  a  number  of  reasons.  The  most 
fundamental  reason  is  that  accurate  ab  initio  predictions  of  the  response  of  compli¬ 
cated  geologic  materials  in  the  regime  of  interest  is  not  now  scientifically  feasible.  If 
it  were,  laboratory  shock  wave  data  would  not  be  as  important  as  it  now  is.  A  major 
reason  for  this  is  that  the  mechanism  for  non-equilibrium  behavior  is  not  well 
understood.  Furthermore,  if  the  theoretical  capability  did  exist,  the  detailed  micro¬ 
scopic  structure  for  the  subject  rock  will  not  be  available.  Finally,  time  constraints 
for  HYDROPLUS  analysis  preclude  detailed  theoretical  analysis. 

The  hysteretic  model,  like  other  aspects  of  the  equation  of  state  is  a  mixture  of 
theory  (thermodynamics)  and  experimental  calibration  (using  shock  wave  data). 

The  shock  wave  data  is  interpreted  in  the  context  of  a  two  phase  model  using 
parametric  models  of  the  low  pressure  phase  (LPP),  the  high  pressure  phase  (HPP) 
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and  the  mixed  phase  Hugoniots.  The  hysteretic  behavior  is  then  partially  derived 
from  thermodynamic  laws,  and  partially  fitted  to  test  data.  In  particular,  the  second 
law  of  thermodynamics  dictates  that  the  fraction  in  the  high  pressure  phase  cannot 
decrease  for  pressures  above  the  LPP-HPP  boundary,  so  the  fraction  in  the  high 
pressure  phase  is  initially  frozen  (held  fixed).  This  procedure  generally  gives 
remarkably  good  agreement  with  the  observed  release  data.  At  lower  pressures 
recovery  of  the  low  pressure  phase  is  possible,  and  the  model  release  curves  are 
adjusted  to  agree  with  observations.  This  fitting  procedure  is  facilitated  by  the 
HYPAWS  EOS  module.  These  issues  will  be  discussed  in  more  detail  in  succeeding 
sections. 
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SECTION  3 

NON-EQUILIBRIUM  CONSIDERATIONS 


Phase  hysteresis  can  be  explained  with  a  simple  physical  model.  The  material 
in  question  begins  in  a  low  pressure  phase  which  is  stable  at  pressures  below  some 
transition  pressure  (Pstart)*  When  the  shock  strength  exceeds  the  transition  pressure, 
the  low  pressure  phase  is  no  longer  stable,  but  the  transformation  to  a  high  pressure 
phase  is  inhibited  by  slow  transformation  kinetics.  The  non-equilibrium  low 
pressxire  phase  deforms  inhomogeneously  and  partially  transforms  to  a  high 
pressure  phase  in  localized  regions.  The  extent  of  transformation  from  the  low 
pressure  to  high  pressure  phases  increases  as  the  shock  pressure  increases.  The 
Hugoniot  state  in  this  transition  region  is  not  in  thermodynamic  equilibrium 
because  of  the  non-equilibrium  phase  composition. 

Consider  an  adiabatic  release  path  from  a  non-equilibrium  Hugoniot  state. 
Non-equilibrium  thermodynamics  constrains  what  can  happen  on  this  path.  As 
long  as  conditions  remain  outside  the  stability  domain  of  the  low  pressure  phase,  no 
reverse  transformation  from  a  high  pressure  phase  to  the  low  pressure  phase  can 
occur.  This  is  a  profoimd  and  inescapable  consequence  of  non-equilibrium 
thermodynamics,  and  leads  directly  to  the  observed  anomalous  releases  below  the 
Hugoniot  in  the  transition  region. 

The  behavior  discussed  above  is  appropriately  called  phase  hysteresis.  The 
inability  to  revert  to  the  low  pressure,  low  density  phase  upon  unloading  has  a 
significant  effect  on  wave  propagation.  The  principal  effect  of  phase  hysteresis  is 
increased  attenuation  of  the  shock  wave.  This  is  because  the  energy  required  to 
transform  from  the  low  pressure  phase  to  the  high  pressure  phase  cannot  be 
recovered  until  the  low  pressure  phase  is  again  stable.  This  is  analogous  to  simple 
models  of  pore  collapse  where  the  energy  involved  in  collapsing  pores  is  never 
recovered. 

The  practical  consequence  of  using  a  hysteretic  phase  transformation  model  is 
that  one  can  easily  calibrate  the  model  to  laboratory  shock  wave  data,  and  then 
accurately  reproduce  pressures,  particle  velocities,  and  times  of  arrival  for  well 
instrumented  underground  nuclear  tests.  Conventional  (i.e.,  non-hysteretic) 
equation  of  state  models  have  failed  to  achieve  this  level  of  agreement  for  pressure, 
particle  velocity,  and  time-of-arrival  simultaneously. 

A  non-equilibrium  equation  of  state  entails  modeling  issues  that  are  not  pre¬ 
sent  for  a  conventional  equation  of  state.  In  general,  the  independent  variables,  mat¬ 
ter  density,  specific  internal  energy  (or  temperature),  and  composition  must  be  speci¬ 
fied  in  order  to  determine  the  pressure.  By  'composition',  we  mean  the  fractions  of 
the  low  and  high  pressure  phases.  For  thermodynamic  equilibrium,  the  composi¬ 
tion  is  imiquely  determined  (from  thermodynamic  theory)  as  a  function  of  the 
matter  density  and  specific  internal  energy.  For  a  non-equilibrium  system,  an 
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independent  criterion  for  determining  the  composition  must  be  provided.  Theory 
provides  very  general  constraints  on  the  composition,  but  very  little  guidance  on 
details.  This  permits  appropriate  empirical  methods  to  satisfy  general  theoretical 
constraints. 

The  detailed  modeling  of  the  non-equilibrium  composition  is  basically 
empirical.  It  should  be  noted  that  this  is  not  essentially  different  than  other  types  of 
equation  of  state  or  constitutive  modeling.  There  is  generally  no  theoretical  basis  for 
choosing  a  particular  functional  form  for  any  constitutive  characteristic.  It  is  up  to 
the  modeler  to  discover  functional  forms  which  satisfy  theoretical  constraints  and 
are  consistent  with  observed  behavior.  Likewise,  the  parameters  of  the  model  are 
chosen  to  agree  with  observations.  The  principal  challenge  to  such  models  is  that 
they  should  continue  to  agree  with  observations  in  new  physical  regimes. 
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SECTION  4 

DISCUSSION  OF  MODEL 


A  computational  non-equilibrium  equation  of  state  which  accounts  for  the  key 
features  of  a  non-equilibrium  silicate  phase  transformation  has  been  implemented 
in  HYPAWS.  This  equation  of  state  is  based  on  a  physical  model  of  the  phase 
transformation  using  the  thermodynamic  mixture  theory  [Gibbs,  1948].  The  model 
assumes  that  the  material  can  exist  in  two  distinct  phases,  a  low  pressure  phase 
(LPP)  which  is  stable  at  low  pressures,  and  a  high  pressure  phase  (HPP)  which  is 
stable  at  high  pressures.  At  intermediate  pressures,  the  two  phases  are  assumed  to 
occur  with  non-equilibrium  concentrations. 

The  HYPAWS  equation  of  state  method  combines  analytic  and  tabular 
components.  The  analytic  part  consists  of  parametric  representations  of  reference 
adiabats  of  the  LPP  and  HPP  constituents,  and  a  model  of  the  mixed  phase  regime. 
The  reference  adiabats  represent  static  loading  from  standard  reference  states  for  the 
low  pressure  and  high  pressure  phases  (room  temperature  and  zero  pressure).  The 
tabular  part  consists  of  a  full  thermodynamic  equation  of  state  which  provides  the 
thermal  pressure  as  a  function  of  the  matter  density  and  thermal  energy.  The 
pressure  is  determined  as  the  sum  of  the  analytic  and  tabular  components.  This  is 
discussed  in  more  detail  below. 

4.1  THE  TABULAR  COMPONENT. 

The  tabular  component  of  the  EOS  consists  of  a  full  thermodynamic  equation 
of  state  for  silicate  rocks  developed  at  S-Cubed.  The  tabular  model  reflects  the 
physics  of  melting,  vaporization,  chemical  dissociation,  and  ionization  and  accounts 
for  the  dependence  of  both  pressure  and  temperature  on  the  material  density  and 
specific  internal  energy.  Due  to  the  complexity  of  the  tabular  component,  it  is  not  ad¬ 
justable  within  HYPAWS.  It  is  also  worth  noting  that  the  tabular  EOS  is  not  just  a 
computational  convenience,  but  is  required  because  the  theoretical  basis  of  the  equa¬ 
tion  of  state  does  not  lead  to  an  analytic  form.  The  following  discussion  provides  a 
brief  summary  of  the  theoretical  basis  for  the  tabular  component  of  the  equation  of 
state. 


The  EOS  is  derived  from  a  Helmholtz  potential  [Callen,  1985],  which 
guarantees  thermodynamic  consistency.  The  Helmholtz  potential  (or  Helmholtz 
free  energy)  depends  explicitly  on  the  temperature  ( ^),  the  density  {P)  and  the 
system  composition  {N-] 

F=F(t?,p,{iVj) 

and  consists  of  a  sum  of  terms  that  account  for  various  physical  components  of  the 
free  energy.  These  terms  are  described  in  detail  by  Kerley. 
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The  formulation  of  the  Helmholtz  potential  is  only  a  first  step  in  constructing 
a  computational  equation  of  state.  At  high  temperatures  where  chemical  dissocia¬ 
tion  or  ionization  occur,  the  composition  of  the  system  is  not  fixed.  If  thermody¬ 
namic  equilibrium  is  assumed,  the  composition,  is  determined  by  minimizing 
the  free  energy  (equivalent  to  maximizing  the  entropy).  Thus,  the  composition  is 
fixed  by  the  temperature  and  density  and  is  no  longer  an  explicit  independent 
variable.  (In  the  non-equilibrium  mixed  phase  regime,  the  phase  composition  must 
be  determined  independently.) 

The  Helmholtz  potential  can  be  decomposed  into  a  cold  contribution  (Fc)  and 
a  thermal  contribution  (Fj) 

F=Fc(p)  +  F^(i?,p) 


The  separation  of  the  computational  equation  of  state  into  an  analytic  term 
that  is  independent  of  temperature,  and  a  tabular  thermal  term  is  related  to  this 
decomposition. 


Given  the  Helmholtz  potential,  thermodynamic  properties  of  computational 
interest  can  be  explicitly  derived,  e.g.,  the  pressure  (P)  and  specific  internal  energy 
(f )  are  given  by 

^  dp 

and 


e  =  F-t?|^ 


Because  P  and  e  are  the  natural  independent  variables  for  hydrodynamic 
computations,  these  equations  must  be  numerically  inverted  to  provide  ^  and  P  as 
a  function  of  P  and  e.  This  presents  a  choice  of  iterating  within  the  calculation  or 
using  tables,  and  is  one  reason  for  using  tabular  equations  of  state.  There  are  other 
reasons  which  are  addressed  below. 


The  Helmholtz  free  energy  applies  to  a  single  phase,  so  the  solid,  liquid  and 
vapor  phases  each  have  distinct  Helmholtz  potentials.  The  equation  of  state  must 
uniquely  accoimt  for  the  phase  composition.  With  the  exception  of  the  solid-solid 
polymorphic  phase  transitions  this  is  accomplished  with  the  classical 
thermodynamic  treatment  of  phase  equilibria.  This  requires  using  the  Gibbs 
potential: 

p 

G  =  F+- 
P 

The  Gibbs  potential  depends  explictly  on  the  pressure  and  temperature.  In 
general,  phase  equilibrium  is  achieved  by  selecting  the  phase  with  the  minimum 
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Gibbs  energy.  Regions  (in  P,  space)  in  which  a  particular  phase  is  stable  are 
separated  by  lines  (coexistence  lines)  for  which  the  Gibbs  energies  in  the  neighboring 
phases  are  equal.  On  the  coexistence  line,  the  matter  density  and  the  specific  internal 
energy  can  vary  continuously  between  those  of  the  bounding  phases.  The  phase 
boundaries  must  be  determined  numerically  with  a  Gibbs/Maxwell  construction 
[Callen,  1985].  This  is  another  reason  for  using  tabular  equations  of  state. 

The  full  equation  of  state  must  accoimt  for  a  number  of  distinct  physical 
phenomena  including:  polymorphic  phase  transitions,  melting,  vaporization, 
dissociation,  and  ionization.  Specialized  first-principles  codes  exist  for  some  of  these 
phenomena  while  others  require  a  more  empirical  approach.  The  complexity  of 
phenomena  and  the  need  to  cover  a  wide  range  of  temperatures  and  pressures 
dictate  that  the  best  physically  based  equation  of  state  models  are  tabular  in  nature. 
Due  to  the  technically  complicated  nature  of  the  high  temperature  equation  of  state, 
it  was  not  feasible  to  include  the  capability  to  adjust  the  thermal  model  in 
HYPAWS.  This  is,  however,  a  potential  area  for  improving  the  HYPAWS 
treatment. 

Shock  propagation  in  the  HYDROPLUS  regime  is  most  sensitive  to  the 
pressure  response  at  relatively  low  temperatures.  Therefore,  it  was  decided  to 
separate  the  full  equation  of  state  into  a  'cold'  component  and  a  thermal  component. 
The  cold  component  is  analytic  and  can  be  adjusted  by  fitting  to  shock  wave  data 
while  the  thermal  component  consists  of  a  fixed  table  derived  from  detailed  high 
temperature  theories.  The  cold  and  thermal  components  are  integrated  into  a 
consistent  model,  and  are  never  applied  in  isolation.  It  is  not  appropriate  to  think  in 
terms  of  a  'thermal  correction'.  The  combined  equation  of  state  depends  on  the 
usual  independent  variables,  matter  density,  specific  internal  energy,  and  in  our  case 
on  the  hysteretic  parameter,  f.  In  addition,  it  also  depends  on  the  parameters  of  the 
analytic  model.  These  parameters  are  adjusted  so  that  the  equation  of  state  agrees 
with  the  shock  wave  experiments.  This  parameter  fitting  is  easily  accomplished 
within  HYPAWS  using  the  EOS  module.  Once  the  parameters  are  set,  the  combined 
equation  of  state  is  used  directly  in  the  finite  difference  solution. 

4.2  THE  ANALYTIC  COMPONENT. 

The  analytic  part  consists  of  parametric  representations  of  adiabatic  loading 
for  the  LPP  and  HPP  using  the  Mumaghan  equation  of  state,  and  a  parametric 
representation  of  the  non-equilibrium  mixed  phase  region.  Three  parameters  are 
required  for  each  pure  phase,  and  four  parameters  for  the  transition  region,  for  a 
total  of  10  parameters.  These  parameters  can  be  set  in  the  HYPAWS  user  interface  to 
obtain  the  best  agreement  between  the  computational  equation  of  state  and 
laboratory  EOS  data.  The  three  parameters  for  the  HPP  and  LPP  are  Po,  Kq,  and  K'o  in 
the  Mumaghan  equation  of  state: 
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Eo 

K 


ff 

p 

-1 

y 

where  is  a  reference  density  (zero  pressure)  for  the  phase  under  consideration,  Kq 
is  a  reference  bulk  modulus  (zero  pressure),  and  K'o  is  the  'slope'  (called  'n'  in  the 
HYPAWS  interface).  From  the  deWtion  of  the  bulk  modulus 


K  =  p 


dp 


the  Murnaghan  equation  gives. 


K  =  K,  +  K'P 


so 


dK 

is  the  slope  of  the  bulk  modulus  as  a  function  of  pressure. 

Generally,  for  rock  K'o  is  in  the  range  3  to  6  and  Kq  is  hundreds  or  thousands 
of  kbar.  The  LPP  and  HPP  phases  are  independently  fit  to  the  full  equation  of  state  to 
determine  the  Murnaghan  parameters. 

The  parameters  of  the  Murnaghan  equation  are  internal  to  the  equation  of 
state  module,  and  are  determined  automatically  when  the  full  equation  of  state  is 
fitted  to  shock  wave  data  using  the  HYPAWS  EOS  module. 

4.3  THE  MIXED  PHASE  REGIME. 

The  solid-solid  pol3anorphic  phase  transition  is  treated  with  an  analytic  non¬ 
equilibrium  model.  The  need  for  the  non-equilibrium  model  is  illustrated  by  con¬ 
sidering  the  laboratory  Hugoniot  (shock)  and  release  data  for  quartz  monzonite 
obtained  by  Ahrens  [personal  communication]  .  Hugoniot  data  are  shown  as  points, 
and  release  data  as  connected  points  in  Figure  4-1.  This  Figure  emphasizes  the  100  to 
500  kbar  range,  where  the  silica  high  pressure  phase  is  stable.  The  four  lower  pres¬ 
sure  release  adiabats  make  a  compelling  case  for  non-equilibrium  mixed  phase 
behavior. 
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Figure  4-1.  Experimental  Hugoniot  and  release  adiabats  for  quartz  monzorute. 
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SECTION  5 

RHYOLITE  ANALYSIS 


Figure  5-1  shows  a  plot  for  SNLA  rhyolite  analogous  to  Figure  4-1  that  was 
available  pretest  for  BEXAR.  This  data  set  is  obviously  sparse  and  problematic.  In 
particular,  there  was  little  data  above  the  quartz  -  stishovite  phase  boundary.  Conse¬ 
quently,  the  equation  of  state  developed  for  the  BEXAR  HYDROPLUS  exercise  had 
to  be  based  on  related  EOS  data  and  judgement  (In  a  real  HYDROPLUS  application, 
much  more  data  is  expected.) 


Rhyolite 


2.5  2.6  2.7  2.8  2.9  3  3.1  3.2  3.3 

Density  (g/cm^ 

Figure  5-1.  Experimental  Hugoniot  and  release  adiabats  for  BEXAR  rhyolite. 

The  densities  obtained  on  shock  loading  show  that  the  silica  is  not  completely 
converted  to  the  high  pressure  phase  until  roughly  400-500  kbar.  (The  LPP  and  HPP 
for  silica  are  quartz  and  stishovite,  respectively.)  Under  static  loa^ng,  this  transfor¬ 
mation  would  be  complete  at  roughly  100  kbar.  This  is  a  clear  demonstration  that 
phase  equilibrium  is  not  established  during  shock  loading.  The  releases  beneath  the 
Hugoniot  are  also  a  consequence  of  the  non-equilbrium  phase  composition,  and 
cannot  be  produced  with  an  equilibrium  equation  of  state. 

To  model  the  non-equilibrium  phase  transition,  the  phase  composition  is  not 
determined  within  the  equation  of  state,  but  by  a  model  (described  below)  that 
determines  the  phase  composition  based  on  the  hydrodynamic  history  of  the 
material. 
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To  deal  with  the  non-equilibrium  phase  transition,  we  must  explicitly 
introduce  the  non-equilibrium  mass  fractions  of  the  low  pressure  phase  (LPP),  fi 
and  the  high  pressure  phase  (HPP),  (2-  The  mass  fractions  satisfy  the  constraint  (mass 
conservation) 

S/.=' 

/ 


It  is  also  required  that  0<fi<l.  The  mass  fraction  of  the  HPP  is  also  represented 
by  f  -  without  a  subscript,  i.e.,  f=f2=l-fi.  f  is  an  additional  dynamical  variable  that 
must  be  tracked  by  the  hydrocode. 

The  material  is  now  modeled  as  a  mixture  of  the  LPP  and  HPP  with  an 
explicit  dependence  on  mass  fraction.  The  specific  volume  (T=l/p)  satisfies 

i 


and  the  specific  internal  energy  (e)satisfies 


where  t,.  and  e,.  are  the  specific  volume  and  specific  internal  energy  of  each  phase  i. 
In  the  Gibbs  representation,  t,-  and  £,-  are  functions  of  P  and  7?,  so  t  and  e  are 
functions  of  P,  7}  and  f.  Note  that  this  is  an  implicit  functional  dependence  because 
P  and  t?  are  functions  of  and  e,..  It  follows  that  P  and  1?  are  functions  of  t,  £  ,  and 
f.  This  is  an  implicit  functional  dependence. 

Once  the  LPP  and  HPP  regimes  have  been  fit,  the  mixed  phase  regime  is 
treated  independently  with  four  additional  parameters.  Two  of  these  parameters 
Pstart  and  Pend/  identify  the  pressures  at  which  the  phase  transition  begins  and  ends 
(the  intersections  of  the  pure  and  mixed  phase  regimes  on  the  Hugoniot)  and  two 
parameters,  nioad  and  nunload/  determine  the  shape  of  the  loading  and  imloading 
curves. 

Pstart/  Pend  and  nioad  explicitly  control  the  shape  of  the  model  Hugoniot  in 
the  mixed  phase  regime.  These  parameters  are  adjusted  so  that  the  model  is 
consistent  with  data  points  that  are  identified  as  mixed  phase  points.  This  is 
illustrated  in  Figure  5-2  with  a  plot  of  shock  velocity  vs.  particle  velocity,  where  the 
mixed  phase  region  extends  roughly  from  1  to  2  km/ s  on  the  abscissa.  All  available 
rhyolite  data  have  been  included  in  this  plot,  and  were  considered  in  developing  the 
model.  The  data  scatter  in  the  low  velocity  region  illustrates  the  need  for  site-specific 
data.  The  fit  in  this  regime  was  not  well  constrained  by  the  data,  and  was  based 
largely  on  scientific  judgement. 
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Figure  5-2.  Shock  velocity  vs.  particle  velocity  for  rhyolite. 

The  incorporation  of  these  parameters  into  the  computational  equation  of 
state  is  transparent  to  the  user,  but  some  further  details  will  be  discussed  for 
completeness.  The  four  mixed  phase  parameters  are  used  in  the  non-equilibrium 
model  for  the  concentration  of  the  HPP,  f.  The  model  for  f  is  constrained  by  the 
second  law  of  thermodynamics,  which  requires  non-negative  changes  in  the 
entropy  as  a  function  of  time.  Therefore,  f  cannot  decrease  unless  P<Pstart/  increase 
unless  P>Pstart-  This  fundamental  requirement  is  responsible  for  the  observed 
phase  hysteresis. 

The  behavior  on  loading  is  thus  given  by 


flood  ~  ®  P^start 


f  toad  ( 


''start 

V  start  ^end  J 


for  Pstart<P<Pend/ 


flood  ~  1  Pend^ 


where  and  are  the  specific  volumes  at  Pstart  Pend/  are  computed  in 
the  code  from  the  respective  fits  to  the  LPP  and  HPP  Hugoniot  data.  Fitting  the 
mixed  phase  Hugoniot  is,  therefore,  equivalent  to  estimating  the  phase  composition 
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on  loading  as  a  function  of  specific  volume.  On  loading,  f=0  or  1  depending  on 
whether  P^start  or  Pend^  pure  LPP  or  HPP). 

On  unloading,  f  cannot  decrease  unless  P<Pstart-  If  we  assume  that  we  have 
converted  completely  to  the  HPP,  we  can  define 


r 


f unload  f 


\^unload 


V  '■ipp 


so  that  f=l  when  t  equals  the  specific  volume  of  the  HPP  at  Pstart/  and  f=0  at  the 
normal  specific  volume  of  the  low  density  phase.  This  expression  assumes  that 
complete  conversion  does  not  occur  without  a  significant  volume  recovery.  The 
parameter  nunload  controls  the  shape  of  this  function. 

Hysteresis,  or  memory,  is  handled  computationally  by  carrying  the 
composition,  f,  as  an  independent  variable  which  can  reflect  the  material  history. 
The  values  of  fioad  and  funload  are  provisional  values  of  f  that  are  determined  by  the 
current  state  alone.  During  loading,  f<fioad/  and  f  is  updated  so  f=fioad/  /  th® 
fraction  in  the  high  pressure  phase  is  increased.  During  unloading  f>fioad/  and  the 
composition  remains  frozen  so  long  as  f<funload-  Finally,  if  f>funload/  f  is  updated  so 
that  f=funload/  i-c-,  the  fraction  in  the  high  pressure  phase  is  reduced.  This  occurs  in 
the  stability  domain  of  the  LPP,  where  the  reverse  transformation  can  occur.  The 
computer  algorithm  for  this  is 

f=MIN(MAX(f,fload),funload) 

Thus,  all  imloads  in  the  mixed  phase  regime  initially  have  a  frozen  phase 
composition  and  then  follow  a  universal  recovery  curve.  There  is  no  independent 
control  of  the  unloads  in  the  frozen  phase  regime.  In  the  final  portion  of  the 
unload,  phase  recovery  is  controlled  by  the  single  parameter,  nunload-  Note  that  it  is 
possible  to  produce  unphysical  results  with  a  poor  choice  of  this  parameter.  It  is  also 
noteworthy  that  this  portion  of  the  unloading  path  is  not  well  constrained  by  data. 
The  value  nunload=l-25  generally  gives  acceptable  results.  It  is  fortunate  that 
uncertainties  in  this  regime  are  not  generally  important  for  ground  shock 
predictions. 

Figure  5-3  illustrates  the  general  behavior  of  the  EOS  model  for  rhyolite.  The 
agreement  between  the  model  and  data  is  illustrated  in  Figure  5-4.  The  monzonite 
model  and  data  are  compared  here  because  of  the  availability  of  data  releases  for 
monzonite  in  the  mixed  phase  regime. 


16 


Pressure  (kbar)  Pressure  (kbar) 


Rhyolite 


Density  (g/cm®) 

Figure  5-3.  Pressure  vs.  density  for  the  rhyolite  EOS  model. 
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Figure  5-4.  Comparison  of  model  EOS  and  data.  The  comparison  here  is  for 
Monzonite. 
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The  parameters  for  the  rhyolite  fit  are: 


Pure  Phase  Parameters 

LPP 

HPP 

Po 

2.41078470  g/cm^ 

4.0  g/cm® 

K„ 

270  kbar 

3.0  Mbar 

n 

4.0 

5.0 

Mixed  Phase  Parameters 

P 

Start 

100  kbar 

Peal 

400  kbar 

n 

load 

1.6 

*^unload 

1.25 

Note  that  the  large  bulk  moduli  lead  to  large  variations  in  pressure  for  small 
variations  in  the  density  making  precision  of  the  reference  density  very  important. 
This  explains  the  high  precision  of  the  reference  density  for  the  LPP. 

We  want  to  emphasize  that  the  hysteretic  unload  is  influenced  by  a  single 
parameter,  nunload/  while  all  other  parameters  associated  with  the  mixed  phase 
region  are  obtained  from  the  Hugoniot  response.  Furthermore,  nunload  has  no  effect 
on  initial  unloading  from  the  Hugoniot  state,  but  is  only  important  at  pressures 
below  Pstart-  It  follows  that  this  unloading  parameter  has  almost  no  effect  on  shock 
evolution  and  yield  estimation.  The  essential  characteristics  of  the  hysteretic  imload 
are  not  modeled  independently,  but  are  physical  consequences  of  the  theory  of  non¬ 
equilibrium  phase  transitions  and  the  parametric  represention  of  the  phase  transi¬ 
tion  on  the  Hugoniot.  The  EOS  prediction  for  initial  release  behavior  near  the 
Hugoniot  thus  represents  an  independent  prediction  of  the  model,  and  is  not 
obtained  by  fitting  release  behavior.  The  comparison  between  the  model  releases 
and  release  data  obtained  is  a  critical  test  of  the  model,  and  the  striking  agreement 
for  monzonite  provides  strong  support  for  this  modeling  approach. 
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SECTION  6 

CONCLUSIONS  AND  RECOMMENDATIONS 


The  development  of  a  practical  methodology  for  treaty  verification  has 
necessitated  an  improved  quantitative  treatment  of  the  dynamic  high  pressure 
behavior  of  geologic  materials.  Non-equilibrium  phase  transformations  in  geologic 
material  have  moved  from  curiosity  and  conjecture  to  clear  fact,  and  theoretically 
sound  models  have  been  developed  to  describe  this  behavior.  Nevertheless,  our 
understanding  of  this  behavior  remains  imperfect.  The  model  presented  here  is 
believed  to  be  adequate  for  the  intended  purpose  inasmuch  as  it  appears  capable  of 
reproducing  shock  loading  and  unloading  data. 

The  model  is  nevertheless  incorrect  in  detail.  In  particular,  investigations 
have  shown  that  the  inhomogeneous  shock  loading  results  in  heated  shear  bands  in 
which  the  transformation  to  the  dense  phase  occurs.  The  heating  clearly  accelerates 
the  transformation  to  the  dense  phase,  and  also  results  in  the  production  of  a  liquid 
phase.  It  is  undoubtedly  an  oversimplification  to  regard  the  high  pressure  phase  as 
pure  stishovite.  Furthermore,  the  possible  role  of  material  strength  has  not  been 
fully  explored. 

In  shelving  HYDROPLUS,  we  are  shelving  an  unfinished  EOS  research 
program.  Improved  models  of  the  shock  induced  phase  transformation  may  emerge 
while  HYDROPLUS  is  on  the  shelf.  Or,  reduced  support  for  science  in  general  may 
result  in  little  advancement  in  this  area.  Should  the  need  for  HYDROPLUS 
reemerge,  there  are  interesting  equation  of  state  issues  which  should  be  reopened. 


r 
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